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Abstract. This paper presents a proposal of a rather new type of effective 
interatomic interaction for molecular dynamics and similar applications. The model 
consists of atoms with prescribed geometric arrangement of active orbitals, represented 
by arms of the length of half of the relevant bond lengths. The interactions have a 
repulsive part between the atomic centers and an attraction between the arm ends of 
different atoms. For each atom pair only the closest pair of arms interacts attempting 
to form the bond. This is a picture of sigma bonds, the pi bonds are modeled by 
alignment of additional internal vectors which might characterize a given atom. Also 
these are primarily pair interactions. Thus there are only pair interactions of several 
types present. The intrinsic arrangement of the model elements, arms and internal 
vectors can be switched depending on the environment. Thus the complexity of the 
many-body potentials is replaced by pair interactions between atoms with complex 
internal behaviour. The proposed model thus allows formation of new geometries, 
establishing new and breaking existing bonds with the use of only pair interactions 
and environment scanning. We discuss in some detail the carbon case and shortly also 
hydrogen, silicon and sulfur. 



Orbital Based Molecular Dynamics 



2 



1. Introduction 

This paper is concerned with a design of a rather new type of effective interatomic 
interaction for molecular dynamics and similar applications. A large body of empirical 
data is concerned with so called bonded force fields, i.e. models where the structures 
in question are unbreakable but deformable, with torsion angles in addition to the 
stretching and compression of the bonds. The model proposed here is of the non- 
bonded type of interactions, also called reactive force fields, which allow formation 
of new geometries and establishing new and breaking existing bonds. During the last 
about 25 years several empirical methods were introduced to represent basic interactions 
between atoms, starting with the three-body interactions of Stillinger and Weber pQ, 
several versions of the Tersoff-Brenner potentials [2], [3], the so called Environment- 
dependent interaction potentials (EDIP) [1], [5], [6] , and the most elaborate ReaxFF 
[7j. These model interactions can be used in molecular dynamics without performing any 
quantum mechanical calculations, they are empirical representations in some cases built 
on the quantum chemical studies, sometimes on models which are fitted to experimental 
facts. The term non-bonded models could also cover the so called ah initio methods 
based on quantum mechanical approach to electronic structure, but these are generally 
considered separately, because their requirments and applications differ from the much 
simpler to implement classical potentials. 

A large group of such non-bonded interactions have been extensively discussed in our 
paper [8], which in fact could serve as a first part of this study. This paper describes a 
proposal for a rather new approach to modeling the interactions between atoms which 
is not based on a single potential function but attempts to model explicitly the known 
geometrical features of bonds and orbitals. In simple terms it can be described as orbital 
based interatomic interaction, which means that the dependence of the energy surface 
on the geometry of atomic configurations is modeled with the help of predefined orbitals. 
We thus wish to call the proposed model OBMD, orbital-based molecular dynamics. It 
will become clear that the model leads to only pair interactions. We refer also to a recent 
work of Rechtsman, Stillinger and Torquato [9 J discussing what they called "synthetic 
diamond and wurtzite structures" which are self-assembled using only isotropic pair 
interactions as an example of special type of interactions. 

Several of the existing potential models already contain features which are beyond simple 
potential functions, as e.g. the environmental dependent interaction potentials (EDIP) 
or the very detailed reactive force fields (ReaxFF). This work has been inspired also 
by these approaches and if the presented framework would become accepted, many 
empirical data could be imported from these models. 

The re-orientation of the orbitals in the presented model in the MD applications would 
lead to the necessity to introduce some model moments of inertia and then the presence 
of spurious kinetic energy associated with the rotation of the atoms with respect to their 
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axes, which clearly does not have any physical meaning. Thus one part of the model 
is also a prescription to avoid this problem by introducing an over-damped rotational 
orientation of the orbitals, much like in the Langevin approach, but here the Langevin 
equations are only applied to the spurious energy which is removed by spurious friction. 

In the typical potential-based applications one is using a potential with predefined 
parameters to simulate a physical situation. From a possible disagreement the 
parameters of the potential can be modified, and the aim is to model the chemistry with 
the basis of the chosen potential functions. Contrary to that, the proposed model aims 
to include as much as possible of the known chemistry by giving the atoms prototypes 
or blueprints of the bonds and the basic geometry as close as possible to the chemical 
bonds, in a way similar to the unbreakable bonds of molecular mechanics but with the 
built in reactivity feature, i.e. the atoms can break and form old and new bonds. At 
the same time the nature of the model is such that the atoms appear to interact only 
via two-body interactions with complex features. 

The described model is not yet implemented in a full width. The parameters needed 
can quite directly be derived from existing model potentials, but this should be followed 
by comparisons and readjustments to ab initio calculations. 

We start by a detailed discussion of a simple planar model in order to introduce the idea 
of the OBMD in section [2j The introduction and discussion of the overdamped angular 
motion is presented in section [31 Modeling of the carbon sheets in three dimensions 
relevant for graphene and graphite is the subject of section HI Extensions needed for 
description of carbon structures, including the possible algorithms for implementation 
of the models of the structural changes are discussed in section [51 Possible applications 
of the proposed approach to other elements and compounds are discussed in section |6j 

2. Simple planar model of the sp 2 hybridized carbon 

The interaction is modelled by overlapping orbitals belonging to each atom (or what we 
imagine as hybridized atomic orbitals). Thus each atom looks more like a little molecule 
with one "heavy" atom - i.e. the actual nucleus, and the centers of the orbitals, which 
are located at the half length of the known bond, as a sort of additional quasiatoms. 
Here we think about the graphene sheet modelled by these objects with three arms. 
The orbitals have maximum overlap when their centers coincide - and when their axes 
are aligned, that is the same as when the nuclei are at the distance of the bond length 
and their "orbital centers" coincide. 

The necessary coordinates are (in 2 dimensions) 



(1) 



The positions of the axes are rigid inside the atom, since we put in the model that the 
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Figure 2. Schematic representation of the model. 



ideal bond length is the property of the two neighbours. The deviation of the actual 
structure from the ideal bond is realized by displacing the two centers and possibly the 
disalignment of the arms, the arms and vertices can only be rotated. 

For characterization of N-atom system are thus needed N atomic coordinates 
supplemented by each atom's orientation angle (in the plane model only) 

{X h Y h A), i = l....N (2) 

where we could also use a more general vector notation Rj = (Xj, Y,j) for the position 
vectors. On each atom the three "arms" and vertices are characterized by their positions 
(vectors) 

2n 



x it a = X,- L + a cos ( ipi + a 

c 

~3 J 



Ui, a = Yi + a sin ( ^ + a^- J , a = 1....3 (3) 
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or 



rj jQ = Rj + a a ijQ a = 1....3 (4) 

where a. ita are the unit vectors of the three arms. For each pair of atoms, their nucleus- 
nucleus relative distance is thus 



r 



{Xi - X 3 f + (Yi - Y t f (5) 



and the nine distances between vertices eq. [3] 

Pia,jf3 = \J (Xi,a ~ Xj^f + (y^ a ~ Vjfif = |r i>a - r jt p\ (6) 

which, however, are functions of only the two atomic coordinate sets eq. [2] (including 
the orientation angle). The distances between the vertices of i-th and j-th atom, pi a ,jp 
must be accompanied by the 9 alignment parameters, or direction cosines explicitely 
given by 

cos 9 iaJI3 = \ [{x i>a - (x jt p - Xj) + (y La - Yi) (y jt/3 - Yj)\ (7) 
i.e. in shorthand notation 

COS 6 iad/ 3 = ai t0l ■ 3L jt p (8) 

For each pair of atoms there is thus only two-body potential, strongly repulsive between 
the nuclei and weakly attractive between the vertices, or orbital centers of the two 
atoms. That is, the nucleus-nucleus repulsion takes the form 

Vn (nj) (9) 
and the strongly coordinating nine-term potential is 

Wij (BU, Rj, ipi, ipj) = 2^ ^2 w (P^ii 3 ' cos QicjP) ( 10 ) 

This all-vertex- vortex interaction is assumed only in this section, and it was the original 
idea of the model, which could be written as a potential. 

In the present model described below, we postulate an additional saturation feature, 
i.e. only the pair of nearest vortices interacts, once the nearest neighbour is established, 
the remaining two on each atom are forgotten, they do not interact. It is a much more 
realistic model, only one electron from each atom is forming the covalent bond (the fact 
that there are so called double bonds will be discussed below). 

To assure that the two vertices form a bond of the required bond length, the two 
arms must point against each other, other configurations must lead to higher energy. 
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The simplest implementation (and thus perhaps not sufficiently flexible for all future 
applications) is a form of w(p, cos 6) containing a sum of an only radial term and one 
product term 

w(p,cos9)=g(p) + f(p)t(cos9) (11) 

where both the radial functions g(p), f(p) and the angular function £(cos#) are strongly 
peaked at p = and cos^ = — 1, respectively. 

In the discussed simple planar model we find that the angular alignment can be 
maintained by the natural stretching resulting from the combination of the nuclear 
repulsion and the attraction g(p) of the vertices, i.e. in 2 dimensions it might be enough 
to consider only the first term of eq. [TT1 without any angular term £(cos0). The above 
addressed feature of saturation, i.e. only one pair of vortices forms the bond, remains 
true also in some three dimensional cases. Thus the atom-atom interaction eq. [TU1 which 
had 9 terms reduces to one term only, 

Wij (Rj , Rj = w (p ia< ijP< , cos 9 ia< dP< ) (12) 

where the two indices a< on i-th atom and (3 < on j-th atom identify the pair of arms 
with the shortest distance between them, i.e. for all a, (3 

Pia < ,jP< — Pia,j/3 

Figure [3] demonstrates that only one pair of atoms can be attached. The third atom 




Figure 3. Excluding the third atom: the filled areas represent the repulsion 
areas of the potential of eq. El 

"learns" that the bond is occupied from the over-all repulsion potential of eq. [9j The 
repulsive part and the attraction between the arms combine together to give a potential 
of the same type as all of the known interatomic potentials. For three and more atoms 
the attraction of the orbitals will be outweighed by the steep repulsion between atoms. 
In the planar geometry the alignment could be established simply by the interplay of the 
repulsion and attraction, without considering the angular alignment function £(cos 9) of 
eq. [TT1 At present we leave open the question whether this additional function is always 
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Figure 4. Left: The interactions of six atoms in random arrangement. Right: 
all ways lead to graphene structure: after the first three atoms have attached 
themselves (A,B,C), the next ones can only fill in gaps (D and E). See also 
figure [3] 

necessary, but clearly its presence will make the fitting of the present model generally 
possible to most of, if not all, the known interactions. 



2.1. Example of a radial shape of the interatomic potential 

It is instructive to consider the interaction energy as a function of distance between two 
nuclei with all the other parameters constant, when the "arms" are kept aligned. This 
radial form of the potential should have a typical intermolecular potential shape. This 
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Figure 5. The radial shape of aligned-arms interatomic potential described by 
eq. [13] with parameters discussed in the text. 

radial form is easily obtained, it is given by the sum of the expression for the repulsion 
potential and the attractive potential of the two aligned arms. If the distance between 
the nuclei is denoted r and the length of each arm is a, then the radial form of V(r) is 
given by 




V(r) = V(r) + g(\r-2a\) 



(13) 
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where V(r) is the repulsion term of eq. Oand g(p) is the radial part of the arm attraction 
potential of eq. [TTJ 

It is quite interesting that by substituting two Gaussian shapes for both V(r) and g(p), 
we always obtain a reasonable shape for this fixed orientation interatomic potential of 
our model. For example, using a=0.7 A and the following two functional shapes 

V(r) = V e~ 3x2 

and 

g(p) = -0.2V e~ 8p2 
we obtain the potential shown in figure [51 



3. Orientation of the bonds: overdamped motion model 

The equations to be obtained from Lagrange function. Schematically: 

dV 

where X is a fictious moment of inertia of the atom and A is also fictious coefficient 
of "angular friction". These two parameters would be needed in the model to assure 
a reasonable behaviour of the simulations. A method which avoids this additional 
complexity is to assume that the angular motion is "overdamped", anologous to the 
over-damped Langevin motion. This means that 

-Aip > Tip 

where A must be adjusted (as a free parameter) and instead of the second order equation 
the angles are updated according to 



4. The graphene sheet model in three dimensions 

In this section we discuss the extension of the previous model to 3 dimensions. The 
model is extended by the rotation axis of each atom, perpendicular to the plane where 
■0 is measured. The unit vector is given by the two customary angles 8 and 0. The 
model thus consists of atoms which contain three arms as in previous section and in 
addition a specified direction vector n which is perpendicular to the plane of the arms. 
The variables describing the atom are thus the arms as in previous section and the 
additional direction vector, specified by the two angles (analogy to the Euler angles). 
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For characterization of N-atom system are thus needed N atomic coordinates R 
supplemented by each atom's three orientation angles 9, ip and ip 

(Ri, 9 u <pi,il>i), i = l....N (16) 

where ip as before specifies the rotation of the arms while the angles 9 and ip specify 
the normal to the plane of the arms, as well as the rotational axis for the rotation by 
ip. Note that we adopt a notation where 

Hi = (sin 9i cos ipi, sin 8i sin <pi, cos 9i) 

The interaction between the arms is identical to the previous section, a new element 
added is the simulation of the 7r-bonding by pairwise alignment of the axes of the 
interacting atoms. 

The nature of the 7r-bonding does not require the saturation feature discussed above, 
rather the requirement that the two atoms engaging in this new interaction also already 
interact via one of the arms-pairs. There are now three angular variables instead of one 
and the equations for angular motion become more complex. The generalized gradient 
and torque will update the angles again through first order equations, since the model 
assumption of over-damped motion is also present here. 

Formally, the interaction can now be written as follows. The atom-atom interaction 
analogous to eq. [JU] becomes now 

w ij ( R i> n i> Rj, n i, 4>j) =w (p ia< j /3< , cos9 ia<ij/3< )+g (\nf ■ nf\ - l) (17) 

where the two indices a < on i-th atom and /5< on j-th atom are given by the minimal 
distance between them ( the pair out of the nine with the shortest distance), 

Pia < ,jP< — Pia,j/3 

for all a, (3 and the Qi a< ,jp < is the angle between these two arms - this time in the full 
three-dimensional space. The vectors nf and nf are the components of the vectors n, 
and iij, perpendicular to the line connecting the two atoms, 

• Tij = nf ■ Tij = 
It should be noted that the term 

fKK-^l- 1 ) 

is a new model effective interaction which can be understood as related to a part of the 
dipole-dipole interaction (cf. [TO]). 

n l ■ n 2 - 2n l ■ n 2 
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Figure 7. Schematic drawing of a group of carbon atoms with sp 2 -hybridized 
orbitals forming a part of the planar graphene structure. 

namely the part which would describe the attraction of two small magnets when they 
are oriented antiparallel to each other and perpendicular to their connecting line. We 
discuss this particular feature and this model mechanism of sheet alignment in more 
detail in a short paper [IT] where also results of a small simulation demonstrate that 
this model interaction really leads to planar alignment. We add that the neglected part 
of the dipole-dipole expression above, i.e. for the vectors aligned with their connecting 
line, can be used to position at the desired distances of the different sheets. Thus the 
representation of graphite situation can be based on split and re-shaped dipole-dipole 
type interaction where the transverse term aligns each sheet while the parallel term 
has such radial dependence with a minimum at the known distance between sheets of 
graphite. 

5. The Carbon Story 

A group of carbon atoms will interact in many different ways, depending on their mutual 
relations. Up to now we have discussed the situation known as sp 2 -hybridization of 
atomic orbitals. In this section we will try to extend this to the other known situations. 
The discussed models should describe all known carbon relations to other carbon and 
hydrogen atoms, including the two diamond structures. 
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The free carbon atom in this model looks as shown in figure [H It is not really a free 




Figure 8. Representation of the so called sp-hybridized orbitals of a carbon 
atom 

atom, which should have electron configuration 2s 2 2p 2 , but it reflects the fact that a 
smallest stray electric field would lead to Stark states in some direction, which then 
would become the two a states, represented by the two arms (with centers of future a 
bonds denoted by small circles). The two (perpendicular) 7r states are shown by the two 
transverse lines perpendicular to the a states and to each other. 

The formula for atom-atom interaction in this case is a simple modification of eq. [T7] 

W'j (Ri, rii, tfj h Rj, n j} ifjj) = w (p ia< ,jp < , cos6 ia<Jf3< )+Y^ 9 (|n^ • nf v \ - l) 

where a < and /3< have the same meaning as in eq. [T71 while the sum runs over the 
representatives fJL, v of the two 7r-orbitals on each atom. The emerging carbon strings 
are referred to as carbyne, and there have been some discussions about its forms with 
alternating single and triple bonds in contrast to equivalent double bonds in a chain. 
It should be noted that the presented model can not distinguish the two alternatives 
unless some further features would be added. 

5.2. Graphene and graphite. 

In a situation when there are more than two close neighbours, it will become energetically 
advantageous to have more electron-orbital mixing than just one Stark mixture, which 
can lead to the formation of additional sigma orbitals (bonds) with neighbour atoms. 
This would switch the atom's behaviour from sp hybridization to the sp 2 hybridization 
shown in figure El with the interaction form already described by eq. [TT1 

In figure El is shown how the sp 2 hybridized states lead to formation of the planar carbon 
structures. 

The 7r-orbitals on all atoms prefer to be aligned, this reflects the formation of 7r-bonds 
and in fact also the origin of the 7r-type electronic energy bands in benzene and in 
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particular graphene. The correlation over many atoms of the 7r-type bands has an 
important energetic effect and this could be built in a future version of the proposed 
model. For the time being we only consider the alignment of the 7r-vectors. 

5.3. Diamond-like structures 

When there are more than three neighbours, the energetic advantage leads to sp 3 
hybridization and geometry shown in figure [9j i.e. there are no orbitals left in the 
pure 7r-form, all four orbitals lead to formation of a-bonds. 

o CJ 4 




Figure 9. Representation of the so called sp 3 -hybridized orbitals of a carbon 
atom, the tetrahedral orbital arrangement 

The formula for atom-atom interaction in this case is again close to the form of eq. [17| 

Wij (Ri , iii , , Rj , nj =w(p ia< J(3< , cos 8 ia< ij/3< ) (19) 

where a< and (3 < have the same meaning as in eq. [T7J but in this case the minimal 
distance is chosen from the possible 16 instead of 9 pairs. Note that now there is no 
dipole-type interaction included. 




Figure 10. Three carbon atoms chain, occuring in diamond, lonsdaleite and in 
aliphatic hydrocarbons. Shown is the plane defined by the three atoms, and the 
three normal planes which contain the remaining orbitals. This arrangement 
is a property of the tetrahedral geometry 

Interaction of tetrahedral sp 3 -atoms would lead to the chain formation of aliphatic 
hydrocarbons in combination with sufficient hydrogen - or to the mesh of carbon layers 
of diamond-like and lonsdaleite-like type structures. 
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The basic three-atom group is sketched in figure [TDJ This three-atom group appears 
both in diamond, lonsdaleite and in hydrocarbons. The bond lengths will be different 
in different situations, in our model the arm lengths are set up by each atom depending 
on its "knowledge" of the neighbourhood. This shows the flexibility of the proposed 
approach. One can put in various instructions for changing the atomic parameters 
depending on the environment, i.e. every atom is a sort of automaton with instruction 
set allowing to enter as much chemistry as desirable. 




Figure 11. The most common isomer of the propane molecule, containing the 
chain of figure [TUl 

One of the situations which can be used to adjust the interaction parameters is 
represented by figure [11] where the hydrogen nucleus is the simple vertex, a line with 
atom at the end, and where the center of the sigma orbital is shown by the circle as in 
the other diagrams, i.e. the molecule of propane. 

5.4- Four body correlations in diamond and lonsdaleite. 

In our previous work [8] we discussed how to include the four-body correlations into the 
Tersoff-type and Stillinger- Weber type models. Here we describe another feature added 
to the present model which will make it possible to energetically differentiate diamond 
from lonsdaleite, i.e. build in an effective four body correlation using the interaction 
form of two-body forces of the type introduced in this paper. 

The two different types of planar carbon chains present in diamond and lonsdaleite are 
shown in figure [12] The chain (a) is already shown in figures [10] and [TTj, and it is present 
in both diamond and lonsdaleite, while the chain (b) is present only in lonsdaleite. We 
shall now construct an effective interaction leading to the alignment of structures of 
the type (a). Comparing the two structures, it is easy to see (especially when using a 
three-dimensional molecular modeling set) that the chains of figure [12] are repeated in 
three planes which have in common the line connecting the two atoms, i.e. one pair 
of the "arms". In the case (a) the three unconnected "arms" on each of the neighbor 
atoms point in exactly opposite directions, it means that in the case (a) all four pairs of 
the arm vectors directionally closest (lying along the same line, but in opposite sense) 
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(a) 



Figure 12. The carbon chains in diamond and lonsdaleite. 



to each other are having scalar products equal to —1. In case (b), and in any other 
conformation of carbon or aliphatic hydrocarbons, this is never the case. The desired 
effective interaction can be described by a spin-spin type formula 



where a iQ , are the unit vectors of the " arms" and the summation runs only over the four 
pairs of the total 16 which are directionally closest, i.e. give the smallest contributions. 
This formula will give zero for the diamond arrangement and a positive energy penalty 
for any other orientation of the two model atoms. When additional atoms are added 
to the ends of the discussed "arms", this two-body interaction results in four-body 
correlations. 

5. 5. Environment dependent switching of interactions 

An essential feature of the proposed model is the switching of a given atom's behaviour 
between the sp 2 to sp l or the sp 3 regimes. In the EDIP model this is implemented 
by adding a more or less continuous function. Here we propose a sort quantum jump 
switching behaviour. However, it does not mean that the structure will change by a 
jump. Only the relation of the arms (and in fact their number, only the a-orbitals are 
represented by arms) If an atom can not find and get aligned with two perpendicular 
pi- vectors, it will switch from the sp 1 -regime of figure [HI to the sp 2 regime of figure El 

The conditions for the switching and their implementation are the strength of the 
proposed model, but on the other hand they represent a problem since at present we 
do not yet have a realistic prescription. This must be obtained both from chemical 
observations and from quantum chemical calculations. 

When considering the sp 3 -regime, we should mention the representation of hydrogen 




(20) 
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atoms, which simply are one-armed atoms with a suitable arm length to complement 
the length of the carbon arm so that their sum becomes the length of the hydrogen to 
carbon bond. 

6. The Silicon and Sulfur Story 

The situation for silicon which has four valence n = 3 electrons is quite different from 
the carbon case of n = 2 electrons. Si-atoms would not experience the sp 2 -hybridization 
of atomic orbitals which is so important for carbon, i.e. clearly no 7r-bonding in addition 
to the sp 3 -hybridization of atomic orbitals. On the other hand, groups of Si atoms can 
arrange themselves in many other geometries (at higher temperatures). The proposed 
model could be used in place of the Tersoff potential for situations where the diamond- 
like structure remains the most important feature. For description of these cases the 
present model concentrating mainly on the features seen for n = 2 atomic states should 
be extended to model also further features of higher atomic states. 

Oxygen and sulfur are another pair of similarly related elements with n — 2 and 
n = 3 orbitals, respectively. These elements show even more different behavior than 
carbon and silicon. In particular the complex behaviour of sulfur aggregates presents a 
problem not well studied in the framework of simple models. Extensions of the present 
model to many such new situations do not seem to be hindered by anything else than 
suitable representations of further features, derived from known chemical observations. 
In particular, the known complexity of sulfur behaviour including the Ss rings (cf the 
work of Stillinger and Weber on liquid sulfur [12]) might be well modelled with the help 
of future extensions to the discussed model. 

7. Conclusion 

We have formulated a procedure for a design of a model empirical atomic interaction 
which has the ability to incorporate many more features than the traditional two- or 
three-body potentials. These many feature are directly connected with observed facts, 
since the interactions are pair interactions only. This is made possible by basing the 
model on predefined orbitals, much like the molecular structure models invented by 
Andre Dreiding. A natural name for the model would thus be the Dreiding interactions, 
but unfortunately the Dreiding force field already exists and does not have any of 
the features discussed here. We have thus chosen the name OBMD for the proposed 
approach. 

The model is not yet practically implemented, we hope that we shall attract experts 
on different aspects to contribute to the fitting and modeling of the details. In 
concluding, we show how powerful the idea is by simply considering the possible shapes 
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of the potentials. In the traditional potentials, there are usually present two different 
inverse powers with opposite sign (Lennard- Jones type), or two exponentials of different 
exponent (Morse type, Tersoff etc), in our approach we can combine even two Gaussians 
to obtain a nice shape for the potential as shown in figure [51 This is naturally made 
possible by the fact that the two gaussians are not centered at the same point, but one 
at the atomic center and the other at the orbital center. 

Acknowledgments 

We would like to thank Dr Norbert Liimmen at University of Bergen for very helpful and 
enlightening discussions of his work based on the ReaxFF system, which contributed in 
a very important way to the formulation of the section on carbon. 

References 

[1] F.H. Stillinger and T. A. Weber, Phys. Rev. B 31, 5262 (1985) 
[2] J. Tersoff Phys. Rev. Lett. 56, 632 - 635 (1986) 
[3] D.W. Brenner, Phys. Rev. B 42, 9458 (1990). 

[4] M.Z. Bazant, E. Kaxiras and J. F. Justo Phys. Rev. B 56 (1997) 8542 

[5] J. F. Justo, M.Z. Bazant, E. Kaxiras V. V. Bulatov S. Yip Phys. Rev. B 58, 2539 - 2550 (1998) 
[6] N. A. Marks, Phys Rev B 63 (2000) 035401 

[7] A.C.T. van Duin, S. Dasgupta, F. Lorant, and W.A. Goddard, J. Phys. Chcm. A 105 (2001) 9396 
[8] L. Kocbach and S. Lubbad, 2009, arXiv:0908.1540, Reactive interatomic potentials and their 
geometrical features 

[9] M.C. Rechtsman, F.H. Stillinger and S. Torquato, Phys. Rev. E 75, 031403 (2007) 
[10] L. Kocbach and S. Lubbad, 2009, arXiv:0908.1548, Geometrical simplification of the dipole-dipole 
interaction formula 

[11] L. Kocbach and S. Lubbad, 2009, arXiv:0908.0092, Transverse dipole-dipole effective interaction 
for sheet arrangements 

[12] F. H. Stillinger and T. A. Weber, J. Phys. Chem. 91, 4899 (1987) Molecular dynamics study of 
chemical reactivity in liquid sulfur 



